load "programs/process.ncl"

begin
    ;定义读取参数
    inpath = "data/"
    outpath = "results/"
    stations = (/"56671"/) ;会理站
    lats = (/26.39/)
    year = 2009    ;2009-2010年会理站干旱过程
    MonthStart = 10
    MonthEnd = 16  ;+6月为次年4月  

    ;计算mci
    data = calculate_mci(inpath, outpath, stations, lats, year, MonthStart, MonthEnd)

    ;计算RMSE
    history = asciiread("data/2009.txt", (/MonthEnd-MonthStart+1,2/), "float")
    result = new(dimsizes(stations), "float")
    time = MonthEnd - MonthStart + 1

    do m=0, dimsizes(stations)-1

        index = data(:, m)
        error = 0.0
        do i=0, time-1;
            if ((history(i,0).gt.index(i)).and.(history(i,1).le.index(i))) then 
                error = error + 0
            else
                distance_left = (index(i) - history(i, 0)) ^ 2
                distance_right = (index(i) - history(i, 1)) ^ 2
                distance = (/distance_left, distance_right/)
                error = error + min(distance)
            end if 
        end do 

        result(m) = sqrt(error / time)
    end do

    ;输出该干旱过程中的RMSE
    print("MCI指数的RMSE为"+avg(result))
end